W poprzednich rozdziałach zdefiniowaliśmy funkcje bazowe na przedziale jednowymiarowym. Pokażemy teraz w jaki sposób definicje tą można uogólnić na problemy dwuwymiarowe. Rozważmy najpierw siatkę prostokątną, na której rozpinamy \( N_x \times N_y=16 \times 16 \) elementów, której jeden brzeg jest równoległy do osi \( x \), a drugi brzeg jest równoległy do osi \( y \). Mamy więc \( N_x=16 \) elementów wzdłuż osi \( x \) i \( N_y=16 \) elementów wzdłuż osi \( y \).

Rozważmy problem aproksymacji bitmapy opisany szczegółowo w rozdziale pierwszym.
Definiujemy następnie wektor węzłów wzdłuż osi \( x \)
[0 0 0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 16 16 16]
oraz wektor węzłów wzdłuż osi \( y \)
[0 0 0 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 16 16 16].
Wektory te definiują funkcje bazowe B-spline trzeciego stopnia o ciągłości dwa, co oznacza, że jeśli policzymy projekcje bitmapy przy użyciu takich funkcji B-spline, nasza aproksymacja będzie bardzo smukła, i możemy policzyć pierwszą i drugą pochodną w każdym punkcje bitmapy. To tak jakbyśmy zrobili bardzo smukły plastikowy odlew naszej bitmapy.

Nasz wektor węzłów wzdłuż osi \( x \) definiuje następujące funkcje bazowe \( B^x_{1;3}(x),...,B^x_{19;3}(x) \), których wzór może zostać wyprowadzony zgodnie ze wzorem Cox-de Boora, ale wyprowadzając B-spline trzeciego stopnia, które rozpościerają się na czterech kolejnych przedziałach. W szczególności będzie ich \( N_x+p=N_x+3=16+3=19 \).
Podobnie wyprowadzamy funkcje B-spline dla wektora węzłów wzdłuż osi \( y \), i dostajemy 19 B-spline'ów trzeciego stopnia \( B^y_{1;3}(x),...,B^y_{19;3}(x) \) zgodnie z formułą .
Następnie nasze dwuwymiarowe funkcje bazowe uzyskujemy przemnażając przez siebie wybraną jedną funkcję bazową z pierwszego zbioru, rozpiętego wzdłuż osi \( x \), oraz jedną funkcję bazową z drugiego zbioru, rozpiętego wzdłuż osi \( y \).
Możemy je uporządkować wierszami (równoległymi do osi \( x \)).
\( B^x_{1;3}B^y_{1,3}, B^x_{2;3}B^y_{1,3}, ..., B^x_{N_x-1;3}B^y_{1,3}, B^x_{N_x;3}B^y_{1,3} \),
\( B^x_{1;3}B^y_{2,3}, B^x_{2;3}B^y_{2,3}, ..., B^x_{N_x-1;3}B^y_{2,3}, B^x_{N_x;3}B^y_{2,3} \),
\( \cdots \),
\( B^x_{1;3}B^y_{N_y-1,3}, B^x_{2;3}B^y_{N_y-1,3}, ..., B^x_{N_x-1;3}B^y_{N_y-1,3}, B^x_{N_x;3}B^y_{N_y-1,3} \),
\( B^x_{1;3}B^y_{N_y,3}, B^x_{2;3}B^y_{N_y,3}, ..., B^x_{N_x;3}B^y_{N_y,3}, B^x_{N_x;3}B^y_{N_y,3} \).
W ten sposób uzyskujemy bazę dwuwymiarowych funkcji bazowych \( \{B_{i,j;3}(x,y)\}_{i=1,...,N_x,j=1,...,N_y}=\{B^x_{i;3}(x)B^y_{j;3}\}_{i=1,...,N_x,j=1,...,N_y}=\{B^x_{i;3}(x)B^y_{j;3}\}_{i=1,...,19,j=1,...,19 } \)
rozpiętych na prostokątnej siatce obliczeniowej zbudowanej z 16*16=256 elementów. Są one przedstawione na Rys. 1.

Możemy rozpiąć na nich aproksymację bitmapy
\( u(x,y) = \sum_{i=1,j=1}^{16,16} u_{i,j} B^x_{i}(x) B^y_{j}(y) \)
obliczając współczynniki \( u_{i,j} \) w sposób opisany w rozdziale pierwszym. W szczególności będziemy musieli wygenerować całki z iloczynów par funkcji B-spline
\( \int_{\Omega} {\color{red}{B^x_{i}(x) B^y_{j}(y)}} {\color{blue}{B^x_{k}(x)B^y_{l}(y)}} dx dy \) dla \( i=1,..,19, j=1,...,19 \). W przypadku funkcji bazowych rozpiętych na płaskiej siatce, całki te możemy rozbić na iloczyn dwóch całek z par B-spline'ów w poszczególnych kierunkach.
\( \int_{\Omega_x} {\color{red}{B^x_{i}(x)B^y_{j}(y)}} {\color{blue}{B^x_{k}(x)B^y_{l}(y)}} dx dy =\int_{\Omega_y} {\color{red}{B^x_{i}(x)}} {\color{blue}{B^x_{k}(x)}} dx {\color{red}{ B^y_{j}(y)}} {\color{blue}{B^y_{l}(y)}} dy \)
gdzie \( \Omega \) to nasza siatka prostokątna, \( \Omega_x \) to brzeg siatki prostokątnej wzdłuż osi \( x \), \( \Omega_y \) to brzeg siatki prostokątnej wzdłuż osi \( y \).

Zajmijmy się teraz ciekawym przypadkiem, w którym nasza siatka obliczeniowa nie jest regularnym prostokątem, tylko rozpięta jest na przykład na rozmaitości, tak jak pokazano na Rys. 2.

W przypadku tym mamy funkcje B-spline liniowe (pierwszego stopnia) w kierunku osi \( x \) oraz funkcje B-spline drugiego stopnia w kierunku osi \( y \). Ilustracją może być tutaj wykres przedstawiony na Rys. 2 to \( y(y-1)(1-x) \).
Wówczas definiujemy nasze wektory węzłów w taki sam sposób, definiujemy jednowymiarowe funkcje bazowe w taki sam sposób, i następnie musimy określić transformacje z siatki prostokątnej na siatkę rozpiętą na brzegu rozmaitości.
Wówczas całki, które musimy wygenerować muszą zostać policzone po brzegu rozmaitości \( \Omega_S \) i musimy inaczej nazwać nasze osie \( \hat{x} \) i \( \hat{y} \) rozpięte na brzegu rozmaitości.
\( \int_{\Omega_S} {\color{red}B^{\hat{x}}_{i;3}(\hat{x}) B^{\hat{y}}_{j;3}(\hat{y})} {\color{blue}B^{\hat{x}}_{k;3}(\hat{x})B^{\hat{y}}_{l;3}(\hat{y})} d\hat{x} d\hat{y} \)
i musimy zmienić zmienne w całkach z prostokąta rozpiętego na brzegu rozmaitości na prostokącie
\( \int_{\Omega_S} {\color{red}{B^x_{i}(\hat{x}) B^y_{j}(\hat{y})}} {\color{blue}{B^x_{k}(\hat{x})B^y_{l}(\hat{y})}} d\hat{x} d\hat{y} = \int_{\Omega} {\color{red}{B^x_{i}(x) B^y_{j}(y)}} {\color{blue}{B^x_{k}(x)B^y_{l}(y)}} |Jac Map(x,y)|dx dy \)
i tutaj pojawia się jakobian odwzorowania \( Map(x,y) \), które przerzuca płaski prostokąt na prostokąt rozpięty na brzegu rozmaitości. Zauważmy, że żeby zastosować algorytm szybkiego solwera zmienno-kierunkowego, jakobian tej mapy musi zostać rozdzielony na iloczyn jakobianów będących funkcją \( x \) i funkcją \( y \), czyli \( Jac Map(x,y)=Jac Map^x(x)Jac Map^y(y) \).
W przeciwnym wypadku nie rozdzielimy całek
\( \int_{\Omega_S} {\color{red}B^x_{i}(\hat{x}) B^y_{j}(\hat{y})} {\color{blue}B^x_{k}(\hat{x})B^y_{l}(\hat{y})} d\hat{x} d\hat{y} = \int_{\Omega} {\color{red}B^x_{i}(x) B^y_{j}(y)} {\color{blue}B^x_{k}(x)B^y_{l}(y)} |Jac Map(x,y)|dx dy =\\= \int_{\Omega} {\color{red}B^x_{i}(x) } {\color{blue}B^x_{k}(x)}|Jac Map^x(x)|dx \int_{\Omega} {\color{red}B^y_{j}(y)} {\color{blue}B^y_{l}(y)} |Jac Map^y(y)| dy \).
Nie jest to zawsze możliwe, i wówczas zamiast solwera zmienno-kierunkowego należy stosować solwer wielo-frontalny lub iteracyjny.